## This is R code for "The Demand for Insurance: Incorporating the Severity of
## Losing Office into the Insurance Model of Judicial Independence"
## Robustness check
## comparing models including component measures

rm(list = ls())
## install.packages("pacman") #load all packages
## remotes::install_github("chrisadolph/simcf") since simcf cannot be directly installed
pacman::p_load(
  tidyverse,this.path,texreg,MASS,AER,ivregEX,ivpack,simcf,arm,stargazer,ivmodel,modelsummary
)
## set working directory
setwd(this.path::here(..=1))
## read data
load("Data/dat.RData")

## Preparation -----------------------------------------------------------------
## a. Select the threats variables
dat <- dat %>%
  mutate(
    threat = threat_cumsum,
    threatS = threatS_cumsum,
    aid_crsc_lag1=ifelse(year>1962&is.na(aid_crsc_lag1), 0, aid_crsc_lag1),
    aid_crsio_lag1=ifelse(is.na(aid_crsio_lag1), 0, aid_crsio_lag1),
    totalaid_lag1=aid_crsc_lag1+aid_crsio_lag1,
    ivvar1 = log((ross_oil_value_2000_lag1 + 1) * (neighbormean_punish_percent_cumsum + 1)+1),
    ivvar2 = log((totalaid_lag1) * (neighbormean_punish_percent_cumsum +1)+1)
  ) %>% #Subset the data to ensure the same observations and comparable BIC/AIC    
  dplyr::select(c(LJI, threatS, threat, gdpS,  demdurS,  country_name, 
                  punish_percent_cumsum, e_h_polcon3_lag1, ross_oil_value_2000_lag1,
                  neighbormean_punish_percent_cumsum, totalaid_lag1, ivvar1, ivvar2)) %>%
  na.omit() %>%
  mutate(interact = (punish_percent_cumsum +1)*e_h_polcon3_lag1,
         interactS = interact/(2*sd(interact,na.rm=TRUE)))

## IV models ------------------------------------------------------------------
iv_total1 <- ivreg(LJI ~ threatS + gdpS + demdurS | gdpS +  demdurS + ivvar1,
                   data= dat); summary(iv_total1,diagnostics=TRUE)

iv_total2 <- ivreg(LJI ~ threatS + gdpS + demdurS | gdpS +  demdurS  +  ivvar1 + ivvar2,
                   data=dat); summary(iv_total2,diagnostics=TRUE)    

iv_total1_components <- ivreg(LJI ~ threatS + I(1+ punish_percent_cumsum) + e_h_polcon3_lag1  + 
                                gdpS + demdurS | 
                                gdpS +  demdurS + 
                                neighbormean_punish_percent_cumsum  + ross_oil_value_2000_lag1 + ivvar1,
                              data=dat);summary(iv_total1_components,diagnostics=TRUE)

iv_total2_components <- ivreg(LJI ~ I(1+ punish_percent_cumsum) + e_h_polcon3_lag1  + 
                                gdpS + demdurS | 
                                gdpS +  demdurS + 
                                neighbormean_punish_percent_cumsum  + ross_oil_value_2000_lag1 + totalaid_lag1,
                              data=dat);summary(iv_total2_components,diagnostics=TRUE)

iv_components_only_oil <- ivreg(LJI ~ I(1+ punish_percent_cumsum) + e_h_polcon3_lag1  + 
                                  gdpS + demdurS | 
                                  gdpS +  demdurS + 
                                  neighbormean_punish_percent_cumsum  + ross_oil_value_2000_lag1,
                                data=dat);summary(iv_components_only_oil,diagnostics=TRUE)

iv_components_only_aid <- ivreg(LJI ~ I(1+ punish_percent_cumsum) + e_h_polcon3_lag1  + 
                                  gdpS + demdurS | 
                                  gdpS +  demdurS + 
                                  neighbormean_punish_percent_cumsum  + totalaid_lag1,
                                data=dat);summary(iv_components_only_aid,diagnostics=TRUE)

## Get values for table --------------------------------------------------------
source("Code/IV_gof_fun.R")

## Get the values for iv_total1
V1.robust.se.pval <- robust_se_pval(dat, iv_total1)
V1.robust.se <- V1.robust.se.pval[1]
V1.pvals <- V1.robust.se.pval[2]
tr1 <- ivreg_gof(iv_total1)

## Get the values for iv_total2
V2.robust.se.pval <- robust_se_pval(dat, iv_total2)
V2.robust.se <- V2.robust.se.pval[1]
V2.pvals <- V2.robust.se.pval[2]
tr2 <- ivreg_gof(iv_total2)

## Get the values for iv_total1_components
V3.robust.se.pval <- robust_se_pval(dat, iv_total1_components)
V3.robust.se <- V3.robust.se.pval[1]
V3.pvals <- V3.robust.se.pval[2]
tr3 <- ivreg_gof(iv_total1_components) 

## Get the values for iv_total2_components
V4.robust.se.pval <- robust_se_pval(dat, iv_total2_components)
V4.robust.se <- V4.robust.se.pval[1]
V4.pvals <- V4.robust.se.pval[2]
tr4 <- ivreg_gof(iv_total2_components) 

## Get the values for iv_components_only_oil
V5.robust.se.pval <- robust_se_pval(dat, iv_components_only_oil)
V5.robust.se <- V5.robust.se.pval[1]
V5.pvals <- V5.robust.se.pval[2]
tr5 <- ivreg_gof(iv_components_only_oil) 

## Get the values for iv_components_only_aid
V6.robust.se.pval <- robust_se_pval(dat, iv_components_only_aid)
V6.robust.se <- V6.robust.se.pval[1]
V6.pvals <- V6.robust.se.pval[2]
tr6 <- ivreg_gof(iv_components_only_aid) 

## save the table  ------------------------------------------------------------
texreg(list(tr1,tr2, tr3, tr4, tr5, tr6),
       file = "./Table/appendix_components.tex",
       override.se = list(V1.robust.se, V2.robust.se, V3.robust.se, V4.robust.se, V5.robust.se, V6.robust.se),
       override.pvalues= list(V1.pvals, V2.pvals, V3.pvals, V4.pvals, V5.pvals, V6.pvals),
       caption="\\textbf{Separating the components of demand for insurance.} 
       Models of de facto independence separating the demand for insurance measure. Models 1 and 2 are fit to the 2SLS models in the main article (there Models 3 and 4), that is fit two one and two instruments, respectively. Model 3 includes the components of demand estimated as separate variables (and instruments for these). Model 4 looks not at the product (i.e. does not use demand for insurance), but only the components individually (also instrumented). Models 5 and 6 break out Model 4, each using only one rather than both of the instruments for competition.",
       omit.coef="name",caption.above=TRUE,
       booktabs=TRUE,dcolumn=TRUE,use.packages=FALSE,
       label="table:iv_alt2",
       custom.coef.names=c("Intercept","Demand for insurance","log(GDP/capita)",
                           "Years democratic (logged)","Percent previously punished","Electoral competition"),
       include.aic=FALSE,include.deviance=FALSE,
       include.rmse=FALSE,style="apsr",digits=2,scalebox=0.75, float.pos="t!") 


## Check the aic bic
## modelsummary(iv_total1, vcov= clustered_vcov_ivtotal1,stars=TRUE)  
## modelsummary(iv_total2, vcov= clustered_vcov_ivtotal2,stars=TRUE)  
## modelsummary(iv_total1_components, vcov= clustered_vcov_total1_components,stars=TRUE)
## modelsummary(iv_total2_components, vcov= clustered_vcov_total2_components, stars=TRUE)  









		